www.gusucode.com > 基于粒子滤波的故障检测,使用似然函数作为检测函数 > 基于粒子滤波的故障检测,使用似然函数作为检测函数/code/FDI based on SIR likelihood/ceshi7.m
close all; clear all; n = 1:600;%sample steps stdw = sqrt(10); npar = 500;%particle number N = length(n);%N denotes the steps of samples H=50:10:200;%threshold falarm=zeros(1,16); missalarm=zeros(1,16);%for every threshold,there is a missalarm rario reponding to it A=50;D=300;F=300;%A denotes the number of simulation,F denotes the number of time points where there is fault j=0; %N=500 % Generate the state process and observations for k=1:16 %threshold from 50,60,70,......200 switch k case 1 h=50; case 2 h=60; case 3 h=70; case 4 h=80; case 5 h=90; case 6 h=100; case 7 h=110; case 8 h=120; case 9 h=130; case 10 h=140; case 11 h=150; case 12 h=160; case 13 h=170; case 14 h=180; case 15 h=190; case 16 h=200; end fN=0; missN=0; %for every threshold,missN denotes the number of missing alarm time points for s=1:50 %for every threshold,50 times simulations is hold on x0 = 0.1; %the prior position of x c0 = 1; b0=25; xpath = zeros(1,N); xmean=0;xvariance=0.1; ymean=0;yvariance=1; xnoise=gauss(xmean,xvariance,N); %the noise of station transite ynoise=gauss(ymean,yvariance,N); b=b0; xpath(1) = x0/2 + b*x0/(1+x0^2) +8*cos(0) +xnoise(1) ; for i=2:N if i<301 b=b0; else b=b0*5; end xpath(i) = xpath(i-1)/2 + b*xpath(i-1)/(1+xpath(i-1)^2) +8*cos(1.2*(i-1)) + xnoise(i); end zpath = 1/20*(xpath.^2) + ynoise; % Particle filter with resampling w = ones(npar,1)/npar; xprev = randn(npar, 1); SParMat = zeros(npar, N); SXParMat = zeros(npar, N); sxparpath = zeros(1,N); likelihood=zeros(N,1); D=zeros(N,1); for i=1:N fnumber=0; missnumber=0; xnext = drawpar(xprev, stdw, i); xs = (xnext.^2)/20; w = w.*(1/sqrt(2*pi)*exp(-((zpath(i)-xs).^2)/2)); l=1/sqrt(2*pi)*exp(-((zpath(i)-xs).^2)/2); Li=sum(l)/npar;%i时刻的所有粒子似然函数值均值 likelihood(i)=Li; Di=0; if i<20 for j=1:i Di=Di+(-(log(likelihood(j)))); end else for j=i-20+1:i Di=Di+(-(log(likelihood(j)))); end end D(i)=Di; if i<301 if D(i)>h fnumber=fnumber+1; else fnumber=fnumber; end else if D(i)<h missnumber=missnumber+1; else missnumber=missnumber; end end w = w/sum(w); SParMat(:,i) = w; SXParMat(:,i) = xnext; sxparpath(i) = w'*xnext; [xprev, w] = impResample(xnext, w); end fN=fN+fnumber; missN=missN+missnumber;%50次仿真中系统故障时残差值小于阈值的时间点总个数 end missalarm(k)=missN/(A*F);%计算每个阈值对应的漏报率,故障漏报率A是仿真总次数F是一次仿真中有故障时的时间点总个数 falarm(1,k)=fN/(A*D); end figure; plot(H,missalarm,'g*-');hold on; %plot(H,falarm,'r--');